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** Abstract : The Van der Pol differential quation is solved by averaging method 
** Subjects: Vibration Mechanics , The Differential equations . 



This worksheet demonstrates Maple's capabilities in finding the graphical solution and dealing with the stability of the steady s 

tate solution of Van der Pol 's differential equation . 

All rights reserved. Copying or transmitting of this material without the permission of the authors is not allowed . 



We consider the Van Der Pol differential equation :_ 

Two topics that we will be in discussion are 

- Finding the steady state solution of this equation by averaging method 

- Estimating the stability of solution obtained . 



x"+CO 2 x - [i (1 - he 2 ) jc' = 
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1 .Define the model Of problem I We examine the effect of non-linear system under external force caused by the 
AC generator ( see fig 1. ) 
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(fig I- ) 

The differential equation of this model is given in the form 



X 

Lx"+(x - a)x'-\ — = E cos vt 
c 



(1) 



2/. .2 



... ,., . , . . x"+k (x -l)x'+x = ecosvt 

Alter simplifying we obtain : v J 

2 . Construct the algorithm . 

Generally the Van Der Pol differential equation can be expressed by 



(2) 



APPROVED 

By CO HONG TRAN at 1:38 pm, Jun 24, 2009 



x"+f(x,x')x'+g(x) = F(t) 



In the special case let f(x,x') = jLl(l-?UC ) and g(x)=CO X 
the equation will be rewritten as : X +CO X — \\,(\ — hx )X = 



By using the transform 



X = *fkx;T = cor; \i x = — 

CO 



(4) 



Then 



x -* co ^" = * ( ° 2 



and substitute these relations to (2) it follows 



co 2 X" ! co 2 X 
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X"+X -\i x (l-X 2 )X'=0 



2\„! 



x"+jc-ji (1-jc z K=0 

Thus we begin with : 

To normalize this equation we find the solution which is expressed in the form 

x = acos ( t + y ) 

It is advantageous to write (p = t + y then the solution will be x = acosq* 

From the transform x= acosty , x' = - asinty . 

We have x" = dx'/dt = -a'sinty - acoscp . (p ' 
By substituting to ( 6 ) , it gives 



(5) 
(6) 



(7) 



-a'sincp -acoscpxp'+acoscp -(0,(1 -a 2 cos 2 (p).(-asincp) = , „, 

In the other hand dx/dt = x' = a'costy - asinty .($>' = -asinty (9) 



Obviously we reach to the following system of differential equations 

-a'sincp -acos(pxp'= -|0.(1- 
a'coscp -<2sin(p.(p'=-asin9 



a'sincp -acoscp.(p'= -\i{\-a 2 cos 2 (p).asin(p -acoscp 



(10) 

> A:=Matrix( [ [-sin (phi) , -a*cos (phi) ] , [cos (phi) , -a*sin(phi) ] ] ) ; 
-sin(tj)) -acos((j))" 



A := 



cos((|)) -asin((|)) 



> u:=vector ( [-mu*a* (1- (a*cos (phi) ) A 2) *sin(phi) -a*cos (phi) , -a*sin(phi) ] ) ; 

u := [-\la (I - a 2 cos((|)) 2 ) sin((|)) - a cos((|)), -a sin((|))] 

> S:=simplify (linsolve (A,u) ) ; 

S := [a (-1 + cos((|)) 2 ) (I (-1 + a 1 cos((|)) 2 ), [i sin((|)) cos((|)) - (I a 1 sin((|)) cos((|)) 3 + 1 ] 

> a(tt) :=simplify (S[l]) ; 

a( tt):= a( -1 + cos( ) 2 ) u. ( -1 + a 1 cos( ) 2 ) 

> phi(tt) :=simplify (S [2]) ; 

())(«) := [lsin((|))cos((|))-(ifl 2 sin((|)) cos((|)) 3 + 1 

By using the symbols a(tt) = a'(t) , §(tt) = §'(t) (11) 

Execute the averaging method for a'(t) and (j)'(^ , we have 

> aO (tt) :=normal (int (mu*a* (l-a A 2*cos (phi) A 2) *sin(phi) A 2 , phi=gamma . .gamma+2*Pi) / (2*Pi) ) ; 
aO(ff) := --|la (-4 + <r) 

o 

> phiO (tt) :=int (mu*a* (l-a A 2*cos (phi) A 2) *sin(phi) *cos (phi) ,phi=gamma. .gamma+2*Pi) / (2*Pi) ; 

0O(«):=O 

The expressions of a ' and (p ' are calculated in the forms : 
a ' = |i.< a(ff):=a(-l + cos((|)) 2 ) ^(-1 +<r cos((|)) 2 ) > 



1 ? 

aO(tt) :=--(! a (-4 + a ) 

o 



(12) 



(p' = |i,< §(tt) := (lsin((|))cos((|))-(ia 2 sin((|)) cos((|)) 3 + 1 > _ <t>0(«):-0 

The steady state solution occurs when aO = or aO = 2 

If aO = then x = x' = this is a trivial solution ( equilibrium ) . 

If aO = 2 then x = 2costy and x' = -2sinq> . 

The necessary and sufficient condition for solution stability includes 

a' =da/dt = *¥(a) with W(ao) = and ¥ 'faoj < 

> Psi(a) :=a0 (tt) ; 

1 i 

¥(a) :=--(! a (-4 + a") 

o 

> DaohamcuaPsi (a) : =normal (dif f (Psi (a) , a) ) ; 

1 3 , 
DaohamcuaPsi( a ) := - (I - - ji a' 

1 o 



(13) 



(14) 



ii ii 



> print ("Dao ham cua ham a' (t) la 

1 3 , 
"Dao ham cua ham a'(t) la : ',' " a"(t)= ",-\i--[ia 2 

2 8 

> subs (a=2 , DaohamcuaPsi (a) ) ; 

-v- 



> print ("Gia tri cua a' ' (t) tai a = 2 la 

"Gia tri cua a"(t) tai a = 2 la : a"(2) = ',' -|i 

> subs (a=0, DaohamcuaPsi (a) ) ; 

1 



a ' ' (t) = " , DaohamcuaPsi (a) ) ; 



a' ' (2) = ", subs (a=2, DaohamcuaPsi (a) )) ; 



> print("Gia tri cua a"(t) tai a = la : a''(0) = ", subs (a=0, DaohamcuaPsi (a) )) ; 

"Gia tri cua a"(t) tai a = la : a"(0) = ',' ■= \i 



Thus if ao = , a"(0) = t-|J. then the solution is not stable 



If ao = , a"(0) = -|i then the solution is stable asymtotically. 

Note : By solving the equation ( 12 ) for the vibration amplitude a(t) (" slowly varying coefficients ") then finding the solution 

expression x(t) of Van Der Pol differential equation , we get 

> aO(tt) ; 

1 2 

--\la(-4 + a ) 

o 

> diff_eq:= dif f (a (t) , t) =-mu*a (t) * (-4+ (a (t) A 2) ) /8; 

a i , 

diff-eq := ^a(f ) = - - \i a(f) (-4 + a(0 l ) 



> init_con:=a (0) =ao; 
init_con := a(0 ) = ao 

> biendo:= [dsolve ({diff_eq, init_con} , {a(t)})]; 

a( = 2 j a( = -2 ■ 



biendo :- 



1 - 



e Qo 2 -4) 
ao 2 



> biendo [ 1 ] : =simplif y (biendo [ 1 ] ) ; 



biendo := a(0 = 2 



1 



2 l-t 1 ') 2 h (~^ f ) 

-ao* + e ao - 4 e 



1 - 



e (ao 2 -4) 
aa 2 



aa 



( accepted ) 



> biendo [2 ] : =simplif y (biendo [2 ] ) ; 



biendo 9 := a(?) = — 2 ■ 



1 



2 , (-M) 2 ^ <"^ f ) 

-ao + e ao" - 4 e 



ao 



( eliminated ) 



> x:=biendo [1] *cos (phi) ; 



x := cos((|)) a(?) = 2 



cos((|)) 



2 H 1 ') 2 ,, <"M 

-ao +e ao -4e 



ao 



with (p = r + y • ( 15 ) 

> x:=t->2*cos (t + gamma) /sqrt( (a A 2-a A 2*exp(-mu*t)+4*exp(-mu*t) ) /a A 2) ; 

cos(? + y) 



x := t-> 2 



2 2 C-f-0 „ <"^" 

a - a e +4e 



The graphical solution of Van Der Pol's equation : 

> a:=0.5; 

a := .5 

> y:=t->2*cos(t + gamma) /sqrt ( (a A 2-a A 2*exp (-0 . l*t) +4*exp (-0 . l*t) ) /a A 2) ; 

cos(? + 7) 

J := t -> 2 t7 

2 2 (-.11) (-.If) 

a* — a~ e + 4 e 



> plot(y(t) ,t=0. .4*Pi) ; 



> animate (2 *cos(x*t + gamma) /sqrt ( (a A 2-a A 2*exp (-0 . l*t) +4*exp (-0 . l*t) ) /a A 2) , x=-6*Pi .. 6*Pi, 
t=0. .10) ; 
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Use graphical method to consider the solution stability of Van Der Pol' s differential equation : 

> with(DEtools) :mu:=0.5; 

> DEplot({ (D@@2) (x) (t)+x(t)-mu*(l-x(t) A 2)*D(x) (t)=0},{x(t) },t=- 

10.. 50, [ [x(0)=l,D(x) (0)=1] ] , stepsize=0.05,title= N Nghiem on dinh cua pt Van Der Pol s ) ; 




3 . Conclusion . 



From graphical results , we reach to conclusion that the steady state solution stability of Van Der Pol diffential equation must be 

precise 

and estimating the property of solution obtained is very necessary . 

As presented above , we might also use the normalization to ( 2 ) by determining the non-trivial solution in the form 

x = Mcoskt + Nsinkt + x* 

(16) 
x' = -kMsinkt + kNcoskt + x *' 



Wl 



th k = 1 



dxWdt = x" = -M'sinr - M cost + N'cost - Nsint + jc* " 



Otherwise 



dx/dt = x' = M'cost -Msint + N'sint + Ncost + x* ' 



= -Msint + Ncost + jc* ' 



Substitute x " to ( 6 ) after simplifying it follows : 

M'cost + N'sint = 

- M's'mt + N' cost = [i[l-(M cost + Nsint) 2 ].[-M s'mt + Ncost] 

> A:=Matrix( [ [cost, sint] , [-sint, cost] ] ) ; 

" cost sint' 
A := 

-sint cost 

> f :=vector( [ [0] , [mu*F] ] ) ; 
/:=[[0],[|iF]] 

> V:=linsolve(A, f) ; 

sint [-[iF] - [0] cos? sint [0] + [-(IF] cast 



V:-- 



sint 2 + cost 2 



sint 2 + cost 2 



We rewrite the expressions : M' = 
N' = [ M-F] c-o^ 



(17) 



With 



F = [I- (M cost + Nsint) 2 ].[-M sint + Ncost] 

Use averaging method for ( 17 ) we get : 

M' = \H<-Fsint> = --M(4-M 2 + 2N 2 -3 N) 



N' = \L< Fcost> = --N (-4-5 M 2 +N) 



(18) 



The steady state solution exists when M = and N = ( trivial solution 
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Or M= 2 or M= -2 and N= 

Thus x = 2cost , x = - 2cost 

If M = and N = or N = 4 Then x = 4sint 

The steady state solutions can be formed generally : 

x = 2cost + 4sinf 
x = -2cosf + 4siru 

But these forms are equivalent to x= 2costy [see (14)] 

Next we begin with Krylov - Bogoliubov approximate method for the Van der Pol's differential equation 

x"+(D 2 x-\i(l-hc 2 )x'=0 

with •* \ x > x ) = — (I ~ "JC ) x , and (J, is a small constant . 

The solution will be estimated by x = r(t) costy(t) 

By taking the first order approximate terms ( neglecting the second and third order errors of constant \i ) we can find the 
amplitude r(t) and the global phase function (p(t) from the following system 
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2% 

ll r r 

r\t) = -^-— f (r cos u-r(d sin u) sin udu = — aAr) 

2710) J 2 

(19> 

(p'(0 = C0H / (r cos w,-rco sin u) cos udu = Ja 2 (r) 

27irco J 



If the initial condition r(0) = ro is satisfied then the solution of ( 15 ) will be equivalent to the solution of second order linear 



x"(t) + a 1 (r o )x\t) + a 2 (r o )x(t) = 



with the error of estimation based on the order of 



(20) 



1^ 



The first order approximate is noticeable in the case of periodic vibration , because the equivalent linear differential equation gives us 
the accumulation and dissipation of energy based on vibration period which we might obtain from the given non-linear differential equation 
Therefore it is useful to apply the equivalent second linear differential equation to observe the non-linear resonant phenomenon . 
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